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1.  Introduction 

We  have  been  working  on  a particular  class  of  interior  point  methods  for  solv- 
ing linear  programming  problems  for  several  years.  (See,  e.g.,  [DBRW91].)  Our 
methods  combine  several  search  directions  that  are  readily  computed  at  each  it- 
eration. The  final  step  is  then  calculated  by  computing  the  step  that  solves  the 
original  problem  restricted  to  the  subspace  spanned  by  these  search  directions.  In 
this  paper  we  propose  an  extension  of  these  ideas  to  the  case  of  convex  quadratic 
programming. 

The  linear  programming  (LP)  problem  that  we  consider  is 

min  u 1 X 

u _ (LI) 

subject  to:  Au  < b 

‘Contribution  of  the  National  Institute  of  Standards  and  Technology  and  not  subject  to 
copyright  in  the  United  States.  This  research  was  supported  in  part  by  ONR  Contract  N-0014- 
87-F0053. 

t Applied  and  Computational  Mathematics  Division,  National  Institute  of  Standards  and 
Technology,  Gaithersburg,  MD  20899.  INTERNET:  boggs@cam.nist.gov 

^Applied  and  Computational  Mathematics  Division,  National  Institute  of  Standards  and 
Technology,  Boulder,  CO  80303-3328.  INTERNET:  domich@cam.nist.gov 

^Applied  and  Computational  Mathematics  Division,  National  Institute  of  Standards  and 
Technology,  Boulder,  CO  80303-3328.  INTERNET:  jrogers@cam.nist.gov 

^Applied  and  Computational  Mathematics  Division,  National  Institute  of  Standards  and 
Technology,  Gaithersburg,  MD  20899.  INTERNET:  witzgall@cam.nist.gov 


I 


April  23,  1991 


and  the  quadratic  programming  (QP)  is 

rmn  ^ u Qu  , . 

subject  to:  Au  < b 

where  u,c  € Q £ A £ and  b £ 3?"”.  We  assume  that  Q is  sym- 

metric and  positive  semi-definite,  i.e.,  we  allow  the  possibility  that  some  variables 
in  the  quadratic  program  enter  only  linearly.  In  addition  we  make  the  standard 
assumptions  that  the  problems  have  a finite  optimal  solution  and  that  A has  full 
column  rank.  We  do  not  assume  that  the  problem  has  a full  dimensional  interior, 
nor  do  we  assume  that  the  problem  is  nondegenerate. 

Interior  point  methods  have  been  extended  to  the  quadratic  programming 
problem  by  several  authors.  For  example,  Shanno,  et  al.  [Sha91]  have  been 
working  on  a primal-dual  method,  and  Anstreicher  has  suggested  a log  barrier 
approach  [AdHRT90]  and  a related  dual  version  [Ans90]  that  is  similar  to  our 
approach,  although  his  work  is  primarily  theoretical.  Our  work  is  based  on  the 
method  of  centers  originally  proposed  by  Huard  [Hua67].  We  believe  that  this 
framework  provides  some  advantages  over  the  primal-dual  approach.  In  particular, 
as  will  be  shown  below,  the  Hessian  matrix  that  is  required  involves  the  sum 
of  Q and  a matrix  of  the  form  Al^ A whereas  in  the  primal-dual  framework  the 
Hessian  involves  a matrix  of  the  form  A^ A.  Clearly,  even  if  Q is  sparse,  the 
inverse  will  generally  not  be,  and  thus  this  latter  form  is  much  more  likely  to 
be  dense.  Also,  since  the  method  of  centers  was  originally  proposed  for  general 
convex  programs,  the  directions  we  compute  are  natural  extensions  of  those  for 
the  linear  programming  case  and  can  be  computed  at  a low  relative  cost. 

In  §2  we  present  a brief  description  of  the  multi-dimensional  method  that  forms 
the  basis  both  for  our  linear  programming  work  and  for  the  quadratic  extensions. 
In  §3  we  discuss  the  generation  of  the  directions;  in  §4  we  describe  the  details  of 
our  implementation;  and  in  §5  we  give  a summary  of  our  preliminary  numerical 
results.  These  results  indicate  that  the  proposed  algorithm  can  solve  the  quadratic 
programming  problems  in  approximately  the  same  time  as  that  required  for  the 
linear  problems.  Finally,  in  §6  we  note  some  directions  for  future  research. 

2.  Multidimensional  Methods 

Our  strategy  for  solving  linear  programming  problems  by  interior-point  methods 
is  based  on  the  observation  that  the  formation  and  factorization  of  the  Hessian 
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matrix  consumes  the  majority  of  the  CPU  time.  Thus  our  goal  is  to  use  the  fac- 
torization to  compute  several  independent  directions,  and  to  determine  the  “best” 
possible  combination  based  on  these  directions.  Viewed  in  this  way,  our  strategy 
is  similar  in  spirit  to  the  “predictor-corrector”  methods  of  Mehrotra  [Meh89]. 

Our  multidimensional  algorithm  for  LP  can  be  thought  of  as  optimizing  over 
a low-dimensional  subspace  at  each  major  iteration.  We  formalize  it  as  follows. 


Multidimensional  Algorithm 

1.  Compute  a strictly  feasible  point,  set  i=0. 

2.  At  u\  generate  q independent  directions,  s\  j = 1, . . . ,^. 

3.  Form  and  solve  the  subproblem 


to  obtain  (*. 


mm  c ' ( + XI  0'^'^ 

j=i 


subject  to:  A j u*  + y~l 
i=i 


< b 


4.  Set  :=  a where  a is  the  steplength. 

5.  If  convergence  then 
Stop; 

Else,  set  i :=  i + 1;  Go  to  2. 


Note  that  the  subproblem  in  Step  3 is  equivalent  to 


mm  > 


^ j=i 
9 

subject  to:  ^ <b  — Au\ 

3 = 1 


The  quantities  c'^sU  As^  ^ and  b—  Au'  need  only  be  computed  once  and  thus  this 
low-dimensional  problem  can  be  solved  efficiently.  We  comment  further  on  its 
solution  in  §4.  The  steplength  a in  step  4 is  fixed  at  .99,  a value  in  accordance 
with  much  of  the  work  in  interior-point  methods  for  LP.  (See  §4  for  details.) 
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The  extension  to  QP  of  this  general  procedure  is  obvious 
objective  function  includes  the  quadratic  term,  i.e., 

min^  + ^('^QgC 

subject  to:  Au  < b 

where 

C,  = + {u‘fQs^ 

and  Qg  G is  defined  analogously. 

In  our  work  to  date,  our  algorithm  optimizes  over  a 3-dimensional  subspace, 
and  hence  we  designate  the  method  by  OSD.  The  subproblem  solution  in  Step  3 
has  three  variables  and  m constraints.  In  the  LP  case,  we  can  efficiently  solve 
this  problem  using  a specialized,  revised  dual  simplex  procedure.  In  the  QP  case, 
this  does  not  appear  to  be  a viable  option,  and  the  QP  subproblem  is  solved  by 
a simplified  interior-point  method,  the  details  of  which  are  in  §4. 

3.  Directions 

The  efficacy  of  the  multidimensional  algorithm  depends  critically  on  the  directions 
that  generate  the  subproblem.  In  this  section,  we  derive  our  directions  from  the 
method  of  centers.  To  do  this,  we  first  present  some  notation. 

Define  the  residuals  for  the  constraints  to  be 

rk{u,t)  = bk  - AkU, 

where  Ak  is  the  A:th  row  of  A.  Define  the  residual  for  the  objective  function  to  be 

ro(u,  t)  = t — c u — \u  Qu 

where  f is  a scalar  whose  value  is  determined  as  follows.  Let  be  strictly  feasible 
and  let  D be  the  value  of  the  objective  function  at  vS.  Note  that  > 

0,  = 1, . . . , m and  that  ro(M°,  f)  > 0 for  t > D.  Thus 

m 

n rk{uO°) 

k=0 

is  positive  in  the  interior  of  the  feasible  region  that  corresponds  to  points  where 
the  objective  function  is  less  than  D.  The  center  can  then  be  defined  as  the 
maximum,  of  this  product,  or,  equivalently, 

m 

min  Z/(u,  = min  V' log  r^(u,  (3.1) 

U U * ^ 

A:=0 


The  subproblem 


(2.1) 
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which  defines  L{u,t).  The  method  of  centers  is  to  solve  (3.1)  to  obtain  set  t^ 
equal  to  the  objective  function  value  at  , and  repeat.  It  can  be  shown  that  this 
yields  a sequence  of  points  that  converges  to  an  optimal  solution. 

In  [DBRW91]  we  show  that  if  the  constraint  on  the  objective  function  is  moved 
continuously  then  a trajectory  is  formed  that  approaches  an  optimal  solution. 
More  importantly,  by  a slight  generalization,  it  can  be  shown  that  there  exists 
a trajectory  connecting  every  strictly  feasible  point  to  an  optimal  solution.  The 
vector  field  defined  by  these  directions  is  given  by 

Our  first  direction,  therefore,  is  known  in  the  literature  as  the  dual-affine 
direction  [ARVS9].  After  some  algebraic  simplifications  we  have  that  is 

(3.2) 


where 


H = A^D^A-\- 


Q 


ro{u,t)' 

Z)  is  a diagonal  matrix  with  A:th  diagonal  entry  l/rk{uA),  and  ffi  is  a scalar. 

It  is  likewise  possible  to  compute  the  so-called  recentering  direction,  i.e.,  the 
Newton  direction  for  solving  (3.1)  given  by 


VuuL[u,t) 

This  direction  is  a combination  of  two  directions,  the  first  of  which  is  and  the 
second  of  which  is  our  second  direction,  5^,  given  by 

De  (3.3) 


where  e is  the  vector  of  all  ones,  and  is  a scalar. 

The  third  direction  is  derived  from  considering  the  rank  one  effect  on  the  Hes- 
sian matrix,  H,  of  the  first  constraint,  whose  index  is  denoted  by  k,  encountered 
in  the  direction.  It  is  easily  shown  that  this  change  is  dominated  by  the  vector 
and  that  the  resulting  direction  is 

(3.4) 


Finally,  using  the  notation  of  factorable  functions  [JM86],  we  obtain  our  fourth 
direction,  a third  order  correction  to  namely 


= H-^EA 


k=l 


T 3 


rk{u,t) 


n2 


(3.5) 
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The  details  for  all  of  these  directions  may  be  found  in  [DBRW91].  Our  strategy, 
outlined  below,  is  to  choose  three  of  these  four  directions  at  each  iteration. 

The  only  difference  between  the  formulas  here  for  the  quadratic  programming 
problem  and  those  for  the  linear  programming  problem  is  the  appearance  of  the  Q 
term.  One  can  readily  observe  that  the  sparsity  of  H is  decreased  by  the  addition 
of  Q,  but  not  catastrophically,  as  it  is  in  methods  such  as  those  mentioned  in  §1. 
The  work  per  iteration,  therefore,  is  approximately  the  same  for  the  quadratic 
program  as  for  the  linear  program. 


4.  Implementation  Details 

In  the  results  reported  here,  we  use  the  basic  03D  algorithm  presented  in  §2  with 
the  three  directions,  s'*,  and  either  or  s^,  to  specify  the  subproblem.  The 
selection  of  the  third  direction  is  based  on  the  proximity  to  the  optimal  vertex. 
The  implementation  uses  in  early  iterations,  and  in  the  final  iterations. 
The  “switch-over”  is  performed  when  the  duality  gap  (see  below)  is  less  than 
or  equal  to  10“'*  or  the  residual  on  the  objective  function  (c.f.,  §3)  is  less  than 
or  equal  to  1/m.  The  main  procedure  continues  until  it  satisfies  at  least  one  of 
three  convergence  criteria:  (a)  the  relative  change  in  the  objective  function,  (b) 
the  relative  difference  between  the  primal  and  dual  objective  values  (see,  e.g., 
[ARV89]),  and  (c)  the  steplength. 

Note  that  computation  of  the  dual  variables  is  theoretically  more  complicated 
in  the  QP  case.  Our  preliminary  work,  however,  simply  extends  the  techniques  we 
used  for  the  LP  case,  using  y = D^A{A^D^A-\-QlrQ)~^[c-[-Qu)  to  approximate  a 
dual  feasible  solution,  provided  that  y > 0.  It  can  be  shown  that  this  is  guaranteed 
to  yield  a dual  feasible  point  in  the  limit,  and  this  appears  to  be  working  well  in 
practice. 

The  subproblem  is  solved  by  using  an  interior  point  approach  on  the  three 
dimensional  subspace.  At  each  iteration  of  the  subproblem,  a dual  afffne  direction 
(3.2)  is  computed  and  a step  is  taken  either  a fixed  percentage  of  the  distance 
to  the  boundary,  or  a distance  which  minimizes  the  objective  function  in  that 
direction,  whichever  is  smaller. 

Problem  scaling,  starting  values  and  the  phase  1 procedure  are  exactly  the 
same  as  used  for  our  earlier  LP  subproblem  work.  In  particular,  both  A and 
the  subproblem  constraint  matrix  defined  in  (1.2)  are  scaled.  Our  algorithm  is 
initialized  by  setting  Uq  = 0,  where  0 denotes  the  0-vector  of  the  appropriate 
dimension,  and  then  taking  a single  recentering  step  using  a quadratic  model  in 
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the  steepest  descent  direction,  AJ De.  When  necessary,  an  initial  feasible  solution 
is  obtained  using  a big-M  Phase  1 procedure  (see,  e.g.,  [BJ77]  and  [DBRW91]). 

The  A matrix  is  stored  in  sparse  format  using  the  XMP  experimental  mathe- 
matical programming  data  structures  described  in  [MarSl].  The  Hessian  matrix 
{A^ D^A  + QItq)  is  stored  using  the  data  structures  from  the  Yale  Sparse  Matrix 
Package  SMPAK  [SMP85].  A minimum- local-fill  ordering  procedure  [VS83]  is  in- 
voked only  once  at  the  beginning  of  the  procedure  to  find  a permutation  of  the 
rows  and  columns  that  reduces  fill-in.  The  Hessian  is  then  factored  and  solved 
using  the  Yale  Sparse  Matrix  Package  SMPAK  [SMP85].  Constraints  that  are 
sufficiently  far  from  the  current  point  u,  are  explicitly  removed  from  the  compu- 
tations. 

The  methods  reported  here  are  implemented  in  Fortran  77  and  executed  in 
double  precision  on  an  IBM  RISC  System/6000  Model  520  workstation  running 
at  20MHz  using  the  x// compiler  with  the  -0  option. 


5.  Results 

We  create  a set  of  QP  test  problems  by  augmenting  83  of  the  linear  programming 
problems  (publicly  available  through  Netlib  [Gay85])  with  a quadratic  term  using 
Q = I.  This  allows  an  easy  comparison  of  the  work  required  to  solve  the  QP 
problems  with  that  required  for  the  LP  problems.  Also,  in  the  sense  that  any 
strictly  convex  QP  can  be  transformed  into  one  with  Q = 7,  this  can  be  regarded 
as  a general  formulation. 

The  convergence  conditions  and  the  corresponding  total  CPU  time  required  for 
these  QP  problems  are  nearly  identical  to  the  LP  results  reported  in  [DBRW91]. 
For  the  83  problems,  40  terminated  due  to  the  convergence  criterion  based  on  the 
relative  difference  between  the  primal  and  dual  objective  values.  Of  the  remaining 
43  problems  that  terminated  due  to  the  objective  function  convergence  criterion, 
all  but  3 found  a point  reasonably  close  (but  not  close  enough  for  duality  gap 
convergence)  to  the  dual  objective  to  indicate  that  the  problem  had  solved  cor- 
rectly. Because  the  “true”  objective  function  values  are  not  readily  available  for 
these  QP  problems,  however,  we  did  not  verify  that  the  correct  objective  value 
was  found.  The  total  CPU  time  for  the  83  QP  problems  is  2637  seconds,  only  50 
seconds  more  than  the  time  required  for  our  best  LP  results. 
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6.  Future  Directions 

Our  preliminary  study  demonstrates  the  feasibility  of  extending  our  multidimen- 
sional method,  OSD,  to  the  QP  case.  We  intend  to  explore  these  ideas  further  by 
first  attempting  the  problem  set  with  a general  positive  semi-definite  matrix  Q. 
In  these  tests,  we  will  not  transform  io  Q = I since  we  believe  that  transforming 
the  problem  io  Q = I will  destroy  too  much  of  the  sparsity.  Next,  although  our 
subproblem  solver  has  performed  adequately  in  these  preliminary  tests,  we  think 
that  some  improvements  are  possible.  Finally,  we  will  investigate  the  use  of  this 
procedure  in  a sequential  quadratic  programming  (SQP)  algorithm  for  general 
large-scale  nonlinear  programming  problems. 
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